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Abstract 


This article investigates the activity regimes of a realistic neuron model 
(as a slow-fast system). The authors study this model using the dynam- 
ical systems theory, for example, qualitative theory methods of slow-fast 
systems. The authors obtain the stability conditions of equilibria in leech 
heart interneurons under defined pharmacological conditions and following 
Hodgkin—Huxley formalism. Although in neuronal models, the membrane 
is usually considered capacitance as a fixed parameter, the membrane ca- 
pacitance parameter is assumed as a control parameter to guarantee the 
existence of Hopf bifurcation using the Routh—Hurwitz criteria. The au- 
thors investigate the transition mechanism between the silent phase and 
tonic spiking mode. Furthermore, some simulations are provided using 
XPPAUT software for analytical results. 
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1 Introduction 


Periodic orbits play an important role in nonlinear dynamical models, espe- 
cially in computational neuroscience from the point of view of their topology. 
Now one of the models that pay attention to researchers, is neuron models 
whose topology is associated with specific and regular activities such as tonic 
spiking, bursting, and resting. 

Most neurons show potential membrane oscillations as endogenous or in 
terms of external perturbations. Many kinds of activities on neuronal models 
based on models following the formalism of Hodgkin—Huxley were described 
in terms of the qualitative theory of fast-slow systems [1, 11]. Fast-slow 
systems are a special class of dynamical systems with at least two distinct 
time scales. The solutions of fast-slow systems possess features that are con- 
strained to stay near the slow-motion manifolds, composed of equilibria and 
periodic orbits of the fast subsystem. If both manifolds are transient for 
the solutions of the corresponding neuron model, then it exhibits a burst- 
ing behavior, a repetitive alternation of tonic spiking and quiescent periods. 
Otherwise, the model indicates the tonic spiking activity if there is a stable 
periodic orbit on the tonic spiking manifold, or it shows no oscillations when 
solutions are attracted to a stable equilibrium state on the quiescent manifold 
[12]. 

Classifying transition mechanisms between various regimes is a fundamen- 
tal problem of the theory of dynamical systems and neuroscience. Real neu- 
rons can show various types of firing patterns, such as tonic spiking, bursting 
oscillations, and silent states, that are frequently observed in neuronal elec- 
trophysiological experiments. We investigate the stability of the silent phase 
and the existence of periodic solutions in the Leech heart interneuron model. 

Malashchenko, Shilnikov, and Cymbalyuk [14] described that the leech 

heartbeat is one of the best-studied invertebrate neuronal networks with 
an identified function, a set of identified participating neurons, and well- 
developed biophysically accurate models. It consists of a small number of 
interneurons distributed over several ganglia. Those located in ganglia 3 and 
4 are responsible for generating basic rhythm [2]. Here we focus on the dy- 
namics of a single interneuron from either ganglion 3 or 4. The canonical 
model has proved itself as a powerful tool for predicting phenomena. 
This model has already been studied by Shilnikov and Cymbaluk [19] and 
Wanga et al. [21]. They consider different parameters, for example, yo ‘or 
Fi, 91, as a control parameter with different initial conditions and investigated 
the existence of different codimension-one or two bifurcations. Researchers 
found the parameter regimes in which bursting, tonic spiking, and silence are 
stable. They located several kinds of multi-stability in which more than one 
activity is stable, such as bi- stability in which bursting coexistent with tonic 
spiking, bursting coexistent with silence, and tonic spiking with silence. 

Cymbalyuk et al. [4] have examined the single heart interneuron model. 
Using extracellular recording methods, they demonstrated that oscillator and 
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premotor cardiac interneurons continue to burst when pharmacologically iso- 
lated with bicuculline, albeit in certain preparations, the bursting is not 
strong. They performed a bifurcation analysis of their model using the leak 
conductance g; and reversal potential E; as parameters. The bifurcation 
diagram shows a narrow strip of parameter values, where bursting occurs, 
separating large zones of tonic spiking and silence. Moreover, the authors 
found multi-stability zones, which are regions where the major activities co- 
exist. They propounded how robustness is achieved in this system, how 
critical variables such as period and duty cycle may be controlled, and how 
modulation may affect this control. 

Shilnikov and Cymbalyuk [19] showed how two different codimension- 
one bifurcations of a saddle-node periodic orbit with homoclinic orbits could 
explain a scenario of transitions between tonic spiking and bursting activities 
in neuron models following Hodgkin—Huxley formalism. They also found the 
co-existence of tonic spiking and bursting modes, which are separated by 
a saddle periodic orbit. This mechanism is based on a Lukyanov—Shilnikov 
bifurcation for a saddle-node periodic orbit with noncentral homoclinic orbits. 
They showed that as a bifurcation parameter is varied, the unstable periodic 
orbit separates the basins of attraction of two co-existing modes: the tonic 
spiking periodic orbit and bursting regime. 

Cymbalyuk and Shilinkov [5] presented conditions under which this model 
demonstrates another kind of bi-stability with two coexistent tonic spiking 
modes. They developed a geometrical framework for Pontryagin’s averag- 
ing method of singularly perturbed systems and located periodic orbits and 
studied their bifurcations. Pontryagin’s averaging method gives a clear geo- 
metrical interpretation of this mechanism. They showed that as a bifurcation 
parameter is varied, the unstable periodic orbit separates the basins of attrac- 
tion of two co-existing modes: tonic spiking regime coexisting with another 
spiking regime. 

Furthermore, by synchronization of Wanga et al. [21] in gap-junction, 
coupled neurons with co-existing attractors of spiking and bursting firings are 
investigated as the coupling strength gets increased for this model. Kolomiets 
and Shilnikov [12] considered this model as an example to investigate generic 
mechanisms of transition between distinct patterns of activity in it. They 
considered neuronal bursting as a modular activity consisting of different 
limiting branches corresponding to oscillatory and equilibrium regimes from 
the fast subsystem, based on the Poincaré return map to analyze nonlocal 
bifurcation to understand better and determine what makes bursting and 
spiking attractors fluctuate their shape and stability. 

The lipid bilayer of the cell membrane is a dielectric that creates a ca- 
pacitance by separating charge inside and outside the cell. Such capacitance 
behavior is represented in electrical equivalent circuit models of the cell mem- 
brane with an electrical capacitor element [16]. The membrane-specific capac- 
itance is considered constant and ubiquitous [6]. Moreover, myriad models 
of excitable cells and axons incorporate a constant membrane capacitance 
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[16]. However, the dielectric properties of living tissues exhibit frequency- 
dependent behavior [7]. 

Howell et al. [10] checked out the effect of a frequency-dependent capacitance 
on neuronal excitability. They found that the effect of dielectric dispersion 
on the sinusoidal stimulation thresholds depended on the electrical properties 
of the neuronal membrane. Moreover, dielectric dispersion increased action 
potential conduction velocity. When analyzing the response qualities of brain 
elements at various frequencies, it is necessary to account for the influence of 
frequency-dependent capacitance . Furthermore, in this paper, we consider 
the membrane capacity as a control parameter in the neuronal model. To the 
best of our knowledge, this is the first work that shows how to occur Hopf 
bifurcation in this model when the membrane capacitance parameter varies. 

The concept of Hopf bifurcation, which is known as the Poincaré— 
Andronov—Hopf bifurcation, seems complicated and abstract. However, this 
theory applies to many real problems. Furthermore, the application of Hopf 
bifurcation theory is quite wide. This theory is applied in tests of various 
oscillations caused by wind gusts (which is very important in the field of con- 
struction), LCR oscillations in electrical circuits, transitions among neuronal 
activity regimes in neuronal models, periodic generation of nerve impulses in 
the nervous system, as well as, for instance in epidemiological models that 
describe fluctuations in the number of patients with an infectious disease. 
Through experimental measurements, the researchers found that the mem- 
brane capacity decreases with increasing frequency. 

In many neuronal models, Hopf bifurcation occurs by changing various pa- 
rameters. Using computer software, the Hopf bifurcation can be determined 
in multidimensional dynamical systems that depend on certain parameters. 
The presence of Hopf bifurcation is proved, and the bifurcation parameter 
value is estimated using the Routh—Hurwitz criterion, which determines the 
difference of this research. 

This manuscript considers a model of a pharmacologically isolated heart- 
beat interneuron when the membrane capacitance is the control parameter. 
Then, the transition mechanism between the silent phase and tonic spiking 
mode is investigated. Hopf bifurcation occurs when a periodic solution or 
boundary cycle, which surrounds the equilibrium point, appears or disap- 
pears with a change in the parameter value. This type of Hopf bifurcation 
occurs where a pair of complex conjugate eigenvalues of the Jacobian matrix 
passes through the imaginary axis while all other eigenvalues have negative 
real parts. Furthermore, the existence of this type of Hopf bifurcation leads 
limit cycle losing stability and disappearing when changing a control param- 
eter. First, we find that such a region for the existence of periodic solutions 
and Hopf bifurcation is being investigated. The type of bifurcation can be de- 
termined analytically, but it is very complicated. However, using computers 
(numerical experiments), we determine the type of bifurcation much faster 
and easier. The trajectories of this model were obtained using Dormand— 
Prince solvers, which are related to Runge-Kutta and are adaptive method 
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of solving ODE. The integration and bifurcation analysis was performed uti- 
lizing Mathematica and XPPAUT. These simulations present our theoretical 
results. 

The manuscript is organized as follows: first, we provide some necessary 
preliminaries. Then we introduce the neuronal model and specify the pa- 
rameters and the initial condition. In Section 3, we investigate a general 
geometrical framework for the analysis of (equilibrium points of a system) 
solutions to slow-fast dynamical systems. In the last section, we present some 
simulations to illustrate analytical results. 


2 Preliminaries 


2.1 Routh—Hurwitz criterion for the existence of Hopf 
bifurcation 


In this section, an important criterion, which is known as the Routh—Hurwitz 
criterion, is presented that gives necessary and sufficient conditions for all 
roots of the characteristic equation that real parts of eigenvalues are nega- 
tive. This criterion is stated in the next theorem; see [20, 17]. 


Theorem 1 (Routh—Hurwitz criterion). Consider a polynomial as follows: 
apa® + Ope" + az_on"-* +++ +a,x+ a0, 


when a > 0, the polynomial has all roots with negative real parts if and only 
if all the leading principal minors of the k x k matrix 


Q2k—-1 A2Q2k-2 °°° Ak 


are positive and a, > 0. We assume that a; = 0 if 7 <0 orj >k. 


Let us consider a system 
t= f(z,u), cER", pweR', (1) 


with an equilibrium (9, uo) and f € C™. 

There is another theorem, which shows the existence of Hopf bifurcation in 
system (1). 

Theorem 2. Assume that system (1) has an equilibrium (9; 19) at which 
the following properties are satisfied: 
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(SH1) Dz f.(xo) has a simple pair of pure imaginary eigenvalues and other 
eigenvalues have negative real parts. 


Thus, there exists a smooth curve of equilibria (x(w); 4) with (uo) = 20. 


The eigenvalues A(js); A(is) of Dz fy.o(x(u)) which are imaginary at u = po 
vary smoothly with ys. In addition, if 


(SH2) d(Re(\(H))/du # 0, 


then a simple Hopf bifurcation will occur. 


Proof. See [8]. 


Moreover, researchers can apply n principal sub-determinants to obtain 
a parameter value for Hopf bifurcation Theorem. It should be noted that 
the characteristic polynomial of the Jacobian matrix J(j1) of system (1) is 
denoted as follows: 


P(A; pm) = det(AIn, — J(u) = po(e) + pi()A + +++ + pa(H)r”", 


where every p;(j) is a smooth function of yz, and p,(s) = 1. We consider the 
case po(1) > 0, because there is not any nonnegative real root. Let 


p2(H) pill 


— 
ee 


pon—1() Don-2(H) : . Pa(ut) 


where p;(u) = 0 if i <0 ori > n. Moreover, when po(1) > 0, according to 
the Routh—Hurwitz criterion, the polynomial P(A; ) of A has all roots with 
negative real parts if and only if the following n principal subdeterminants 
of L,() are positive: 


¢ Dj(u) = det(Li(u)) = pi(u) > 9, 


+ Daly) = det(La(u)) = det (PM PU) 5 o, 


© Dy(u) = det(Ln(u)) > 0. 


Since Dp (4) = Pn({t)Dn—1(“) and for characteristic equation pp (i) = 1, the 
Routh—Hurwitz criterion conditions can be formulated as 


Po(L) >0,Di > 0, D2 >0,-.-,Dn—-1 > 0. 
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Theorem 3. Suppose that there is a smooth curve of equilibria (x(), 1) 
with x(49) = 2 for (1). Then conditions (SH1) and (SH2) for a simple 
Hopf bifurcation are equivalent to the following conditions on the coefficients 
of the characteristic polynomial P(); /4): 


1. po(to) > 0; Dy (uo) > 0;...; Dn—2(H0) > 0; Dn—1(H0) = 0, 
2. Poe ttin) #0. 
Proof. see [13]. 


2.2 Neuron model 


Because a complete study of a neuronal model including all the currents iden- 
tified experimentally is very complex, we study the dynamics of a reduction 
of the canonical leech heart interneuron model. Depending on the parame- 
ter values, the model may display several regimes, including tonic spiking, 
bursting, and quiet phase. These regimes represent the intricacy of the dy- 
namics of diverse membrane ionic currents acting on distinct time scales. 
These neurons are particularly attractive for analysis in terms of the theory 
of dynamical systems since their membrane ionic currents were measured via 
voltage-clamp experiments and well described by a canonical model using 
the Hodgkin—Huxley formalism. We focus on the dynamics of a single in- 
terneuron when isolated pharmacologically from the rest of the network. In 
these neurons, eight voltage-dependent ionic currents have been well identi- 
fied and characterized [18, 9]. These currents are separated into four groups 
and classified by their ionic specificity in Table 1. Moreover, “C.N.” is the 
abbreviation of current names. The model equations for Iy, current were 


Table 1: Ionic currents in leech heart interneuron model 


Ton C.N. Description Method of quantify 
Ina fast sodium current Hodgkin—Huxley 

Na é : 
InaP persistent sodium current Voltage clamp 
Ti delayed rectifier-like potassium current 

K Ika persistent potassium current Voltage clamp 
Ika fast transient potassium 

Ca Icar rapidly inactivating current Voltage canis 
Ioas slowly inactivating current 

Na, K In hyperpolarization activated current Voltage clamp 


borrowed from the original work by Hodgkin and Huxley, which adjusted 
for leech kinetics. These currents are independent of the intracellular con- 
centration of any particular ion. A canonical model of a single neuron was 
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described by a system of 14 differential equations running at multiple time 
scales which vary from a few milliseconds through seconds. As pointed out 
before, a complete investigation of this concept would be exceedingly tough 
and difficult. By blocking clusters of currents in live cardiac interneurons and 
exploiting their distinctive characteristics, neuronal dynamics are simplified. 
From the standpoint of the theory of dynamical systems, these distinctive 
behaviors provide intriguing phenomena for investigation [19]. 

Shilnikov and Cymbalyuk [19] removed from the 14-D canonical model de- 
veloped in [9] the equations and terms describing blocked currents: Icar, 
Icas, and I;,. For simplicity, they assumed that the partial block of outward 
currents completely removes Ix, as well as Ix4, and the current Iyap is 
ignored, whereas it reduces Ix2. Now, the neuron can be described by the 
model based on just two currents activation of [x2 as mx and inactivation 
of Ing as hnag and V as the membrane potential. Therefore the resulting 
model described in [3] is as follows: 


CV! = —(Gx2Migg(V — Ex) + gi(V — Ex) 
+9naf (—150, 0.0305, V)?hna(V _ ENa)), 


ee f(—83,0.018+-V 22°F" V)—mxe (2) 
K2—_ TK2 ’ 

n! _ f(500,0.03391,V)—hnva 
Na TNa 5) 


where the variables V, mx2, and hy display as membrane potential, acti- 
vation of Ig, and an inactivation of Iy,, respectively. 


Moreover, f is a Boltzman function: f(#,y,z) = and other 


1 
(pert) 


Table 2: Electrical properties of the leech heart interneuron Model 

Parameter Description Unit(s) | Value 
Cc Membrane capacitance nF B.P. 
Ina Maximum conductance of Inq nS 200 
9K2 Maximum conductance of Ix nS 30 

I Leakage conductance of current ns 8 
Ena Sodium reversal voltage mV 45 
Exe Potassium reversal voltage mV -70 
Ey Leakage reversal voltage mV -46 
Shift of the membrane potential 
setae a of half-inactivation of Ix2 mV -25.98 
from its canonical value 

TNa Time constants of inactivation of In, ) 0.0405 

TK2 Time constants of activation of Ig S 0.9 


parameters are described in Table 2. As we noted before, we use C' as a 
bifurcation parameter or (B.P.) in system (2). 
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3 Main results 


Now, we replace the parameter values are given in Table 2 in system (2) and 
compute its equilibrium points. These points are as follows: 


a = (0.004327, 0.424776,0), gz = (—0.02793, 0.04831, 0.04788). 


To investigate the dynamical behavior of equilibrium points, we need to 
get the Jacobian matrix of system (2) at these points. For this aim, we obtain 
the Jacobian matrix of system (2) as follows: 


90000. —0.045 200(2—0.045 
can cea J2 _ 30y2-—8 —60(a +0.07)y eas ) 
92.2222C'e—83(#—0.00798) “141111¢ 0 


(e—83(@—0.00798) 41 2 ; 


12345.7Ce00(#+0.03391) 
(€500(@+0.03391) +1)? 0 —24.6914C 


where (V,mx2,hNa) = (2, y,2) and a = e~150(+0.0305) | The Jacobian ma- 


trix and the characteristic equation of system 2 at q,; and q2 are computed 
below 


e at qm 
—13.4131 1.89443 8.00404 
A = 22.5345C —1.11111C 0 
—0.0000613486C 0 —24.6914C 
P,, (A, C) = 1422.06C? + (27.4349C? + 388.782C) 
+(13.4131 + 25.8025C)\? + d3, 
e at qe 
16.7353 —0.121944 3.07555 
A,(c) = | 4.24029C —1.11111C 0 
—562.805C 0 —24.6914C 


Py, (A, ¢) = 1476.9C? + (1299.64C + 27.4349C7)\ 
+(25.8025C — 16.7353)\? + A°. 


Now, we investigate the stability of equilibria by the Routh—Hurwitz cri- 
terion. Therefore we need to compute L;(C),i = 1,2 and find where the 
coefficients of the characteristic polynomials are positive. In the following, 
we denote the coefficients of P,,(A,C) by Pi;, where j shows the power of 
parameter. Hence the Routh—Hurwitz criterion holds for P,, (A, C) if 
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Pio = 1422.06C? > 0; 

Py, = 27.4349C? + 388.782C > 0; 
Pry = 25.8025C + 13.4654 > 0; 
Pi3=1>0; 


* D,(C) = det(Ly(C)) = pi(C) = 27.4349C? + 388.782C > 0; 


* D(C) = det(L2(C)) 
— det ( 27-4349C? + 388.7820 1422.06C2 af 
=a 1 25.8025C + 13.4654 


We have similar conditions for P,, (A, C) as follows: 
Poo = 1476.9C? > 0; 
Po, = 27.4349C? + 1299.64C > 0; 


P22 = 25.8025C — 16.7353 > 0; 
Po3 = 1>0; 


¢ D,(C) = det(L,(C)) = pi(C) = 27.434907 + 1299.64 > 0; 


* Do(C) = det(L2(C)) 
— det ( 27-4349? + 1299.64C 1476.90? = 
=e 1 25.8025C — 16.7353 


By the above computations, we have the following theorem. 


Theorem 4. Let q, = (0.00433, 0.42478, 0) and gz = (—0.02793, 0.04831, 0.04788) 
be the equilibrium points of system (2). Then 


1. q is unstable for 0 < C <1, 


2. gg is stable for C > 0.678033 and unstable for C < 0.678033. 


Proof. It is trivial. 


If we put Co = 0.678033 at gq. = (—0.02793, 0.04831, 0.04788), then we 
get the following relations: 


Po (0.678033) = 1476.9(0.6780332) > 0; 
P>; (0.678033) = 27.4349(0.6780332) + 1299.64(0.678033) > 0; 


27.4349CZ + 1299.64C 1476.9C%4 
1 25.8025Co — 16.7353 


= 707.889C% + 31597.9Cp — 21749.9 ~ 0; 


D2(0.678033) = det ( 
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dD2(C) d(707.889C? + 31597.9C' — 21749.9 


) 
dC C=Co dC cig iranag 
= 1415.778(0.678033) + 31597.9 4 0. 


Therefore, according to these results we can state the following theorem. 


Theorem 5. Suppose that for Co = 678033 and equilibrium point g = q2 of 
system (2) the following relations satisfy: 


2. Po; = 27.4349C0? + 1299.64Cy > 0; 


dD2(C) 


Then Co is a simple Hopf bifurcation value for system (2) at the equilibrium 
point q = q. 


4 Numerical simulation 


In this section, several numerical simulations are presented to verify the effi- 
ciency of our analytical results. In addition, we draw the bifurcation diagram 
that shows periodic orbits crossing the plane V = —0.025. The simulation 
parameter values are shown in Table 3. We study the stability of the system’s 
equilibrium points for particular C parameter values satisfying Theorems 4 
and 5. 

If we consider the starting point p,; = (0.0043, 0.42,0.0001) near equilib- 


Table 3: C parameters values and initial values 


C Initial value Figure 
0.64 (0.004,0.413,0.001) Figure 1 
0.685 | (-0.0266,0.0478,0.0475) | Figure 2 

0.6 | (-0.0266,0.0478,0.0475) | Figure 3 


rium point qi, for C = 0.64, then it is observed that the equilibrium point is 
unstable; see Figure 1. 
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(a) (b) 
: \\ 
(c) (d) 


Figure 1: (a). Phase portrait in the space (V,mxK2,hna) at (0.0043, 0.42,0.0001) and 
C = 0.64. (b), (c), and (d). Time history of the neuronal model in terms of V, mxK2 and 
hna, respectively. 


(a) (b) 


(c) (d) 


Figure 2: (a). Phase portrait in the space (V,mxK2,hNna) at (—0.0266, 0.0478, 0.0475) and 
C = 0.685. (b), (c), and (d). Time history of the neuronal model in terms of V,mxK2 
and hyna, respectively. 
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Now suppose that the starting point po = (—0.0266, 0.0478, 0.0475) near 
equilibrium point go, for parameter C' = 0.64. Then it is observed that the 
equilibrium point is unstable; see Figure 2. 

Therefore, by Theorem 5, Co = 0.678033 is a Hopf bifurcation value 


Pili wo 
i" is 
h f 

“| 


(c) (d) 


Figure 3: (a). Phase portrait in the space (V, mK2,hna) at (—0.0266, 0.0478, 0.0475) and 
C = 0.6. (b), (c), and (d). Oscillatory waveforms generated by the neuron model (2) in 
terms of V,mxK2 and hna, respectively. 


for equilibrium gz of leech heart interneuron model. Furthermore, for 
Co, the Jacobian matrix A,,(Co) has a pair of pure imaginary eigenvalues 
Ay; A2 & +29.8422097, and a negative real eigenvalues 43 = —0.762921. 

We solved equation (2) by means of the fourth and fifth-order Runge-Kutta 
(RK-45) method. To show the occurrence of Hopf bifurcation, we present 
some numerical simulations. The related trajectories in the phase space 
(V, m2, hn) and the oscillatory waveforms generated by the neuronal model 
(2) for different values of parameter C' are presented in Figures 3 and 4. 

In Figure 5, we show that by considering a fixed point and changing C, the 
periodic orbit at the bifurcation value disappeared. As it is illustrated in 
Figure 4, the equilibrium point 


g2 = (—0.02793, 0.04831, 0.04788) 
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HII 


arnilll| II| \| 


“= Mb) 


Figure 4: Comparison between phase portraits and oscillatory waveforms in (V,hwa)- 
space for C = 0.645 and C = 0.72. 


mk 


Figure 5: Bifurcation diagram of system (2) 


is a stable focus when the membrane capacitance is bigger than Co. In this 
case, through time, the amplitude of oscillatory waveforms will vanish. While 
for C < Co, the equilibrium point g2 turns out to be an unstable focus 
surrounded by a stable limit cycle. 


This subsystem exhibits either tonic spiking, which corresponds to a stable 
limit cycle or quiescence which corresponds to a stable equilibrium point. 
Stability loss of the stable limit cycle in the fast subsystem gives rise to a 
stable equilibrium through a supercritical Andronov—Hopf bifurcation. When 
a stable limit cycle shrinks into an unstable equilibrium, then the equilibrium 
becomes stable. That is, the equilibrium undergoes supercritical Andronov— 
Hopf bifurcation. 
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5 Conclusion 


Since finding bifurcations, such as Hopf bifurcations in dynamical systems, is 
not an easy task, mathematical researchers are looking for computer tools and 
methods to observe Hopf bifurcations and estimate the number of bifurcation 
parameters. 

In this paper, a leech heart interneurons model has been considered when 
membrane capacitance is the control parameter. Then by using the Routh— 
Hurwitz criteria, the sufficient conditions have been examined such that under 
which the system undergoes supercritical Andronov—Hopf bifurcation. Hence, 
in this model, the transition between the silent phase and tonic spiking mode 
was through Hopf bifurcation. Moreover, the authors introduced an interval 
in the parameters space in which the system has a stable limit cycle in which 
the neuron is in a tonic spiking mode. Numerical simulations were carried 
out using the fourth and fifth-order Runge-Kutta method to support our 
analytical results. In fact, by simulations and Hopf bifurcation theory, we 
showed that the stable limit cycle could be disappeared. 


Acknowledgements 


Authors are grateful to there anonymous referees and editor for their con- 
structive comments. 


References 


(1] Bertram, R., Butte, M.J., Kiemel, T. and Sherman, A. Topological and 
phenomenologial classification of bursting oscillations, Bull. Math. Biol. 
57(3) (1995) 413-439. 


(2) Calabrese, R.L., Nadim, F. and Olsen, @.H. Heartbeat control in the 
medicinal leech: a model system for understanding the origin, coordi- 
nation, and modulation of rhythmic motor patterns, J. Neurobiol. 27(3) 
(1995) 390-402. 


[3] Cymbalyuk, G.S. and Calabrese, R.L. A model of slow plateau-like oscil- 
lations based upon the fast Na+ current in a window mode, Neurocom- 
puting 38 (2001) 159-166. 


[4] Cymbalyuk, G.S., Gaudry, Q., Masino, M.A., and Calabrese, R.J. Burst- 
ing in leech heart interneurons: cell-autonomous and network-based 
mechanisms, J. Neurosci. 22 (24) (2002) 10580-10592. Neurosci. 22, 
10580 (2002). 


Iran. j. numer. anal. optim., Vol. 13, No. 2, 2023,pp 170-186 


185 


Stability and Hopf Bifurcation in Leech Heart Interneuron Model 


[5] Cymbalyuk, G. and Shilnikov, A. Coexistence of tonic spiking oscilla- 


13 


14 


15 


[16 


[17 


[18 


tions in a leech neuron model, J. Comput. Neurosci. 18(3) (2005) 255— 
263. 


Gentet, L.J., Stuart G.J., and Clements J.D. Direct measurement of 
specific membrane capacitance in neurons, Biophys. J. 79 (2000) 314— 
320. 


Grimnes, S. and Martinsen, @.G. Alpha-dispersion in human tissue, J. 
Phys.: Conf. Ser. 224(1) (2010) 012073. 


Guckenheimer, J. and Holmes, P.J. Nonlinear oscillations, dynamical 
systems and bifurcations of vector fields, Vol. 42. Springer Science & 
Business Media, 2013. 


Hill, J., Lu, M., Masino, O. Olsen, R. and Calabrese, L. A model of a 
segmental oscillator in the leech heartbeat neuronal network, J. Comput. 
Neuroscience. 10 (2001) 281-302. 


Howell, B., Medina, L.E. and Grill, W.M. Effects of frequency-dependent 
membrane capacitance on neural excitability, J. Neural Eng. 12 (2015) 
056015. 


Izhikevich, E. Neural excitability, spiking and bursting, Int. J. Bifurc. 
Chaos 10 (6) (2000) 1171-1266. 


Kolomiets, M. L. and Shilnikov, A. L. Poincaré return maps in neu- 
ral dynamics: three examples, International Conference on Difference 
Equations and Applications, pp. 45-57. Springer, Cham, 2019. 


Liu, W. Criterion of Hopf bifurcation without using eigenvalues, J. Math. 
Anal. Appl. 182(1) (1994) 250-256. 


Malashchenko, T., Shilnikov, A. and Cymbalyuk, G. Bistability of burst- 
ing and silence regimes in a model of a leech heart interneuron, Phys. 
Rev. E 84(4) (2011) 041910. 


McIntyre, C.C., Richardson A.G., and Grill, W.M. Modeling the ez- 
citability of mammalian nerve fibers: influence of after potentials on the 
recovery cycle J. Neurophysiol. 87 (2002) 995-1006. 


McNeal, D. R. Analysis of a model for excitation of myelinated nerve 
IEEE Trans. Biomed. Eng. 23 (1976) 329-337 . 


Monfared, Z. and Dadi, Z. Analysing panel flutter in supersonic flow by 
Hopf bifurcation, Iranian Journal of Numerical Analysis and Optimiza- 
tion 4(2) (2014) 1-14. 


Opdyke, C.A. and Calabrese, R.L. A persistent sodium current con- 
tributes to oscillatory activity in heart interneurons of the medicinal 


leech, J. Comp. Physiol. 175. (1994) 781-789. 


Iran. j. numer. anal. optim., Vol. 13, No. 2, 2023,pp 170-186 


Parvizi, Razvan and Alipour Fakhri 186 


[19] Shilnikov, A.L. and Cymbalyuk G. Transition between tonic spiking and 
bursting in a neuron model via the blue-sky catastrophe, Phys. Rev. Lett. 
94 (4) (2005) 048101. 


[20] Siili E. Numerical solution of ordinary differential equations, Mathemat- 
ical Institute, University of Oxford, 2010. 


[21] Wanga Q., Duana, Z., Fengc, Z., Chena, G. and Lu, Q. Synchronization 
transition in gap-junction-coupled leech neurons, Phys. A: Stat. Mech. 
Appl. 387 (2008) 4404-4410. 


How to cite this article 

Parvizi, F., Razvan, M. and Alipour Fakhri, Y., Stability and Hopf Bifurca- 
tion in Leech Heart Interneuron Model. Iran. j. numer. anal. optim., 2023; 
13(2): 170-186. https: //doi.org/10.22067 /ijnao.2022.75087.1101 


Iran. j. numer. anal. optim., Vol. 13, No. 2, 2023,pp 170-186 


